%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This Matlab program is for 'The explicit formula for the Hodrick-Prescott filter in finite sample' by Adriana Cornea-Madeira
%This probram compares the results from Theorem 1 with the HP filter based on the direct inversion of the 
%matrix in (2)
clear all; clc;

n = 12; % sample size
dates = (1:1:n)';
c = randn(n,1); 
tau = zeros(n,1); eta = randn(n,1);
tau(1) = eta(1);
for i = 2:n
    tau(i) = tau(i-1) + eta(i);
end

%data
Y = tau + c; 

% smoothing parameter
alpha = 5; 

i = linspace(1,n,n)';
j = linspace(1,n,n)';
T = sin(i*j'*pi/(n+1));
la = 1+alpha*(2-2*cos(i*pi/(n+1))).^2;
La = diag(la);
Lai = diag(1./la);
T = sqrt(2)*T/sqrt(n+1); %matrix of eigenvectors
S = sqrt((n+1)/2);

% computation of matrix K1 with typical element (7)
i = 1:2:n;
j = 1:2:n;
K1 = zeros(n,n);
for k = 1:length(i)
    for m = 1:length(j)
        K1(i(k),j(m)) = (2*sin(i(k)*pi/(n+1))/S-sin(2*i(k)*pi/(n+1))/S)*(2*sin(j(m)*pi/(n+1))/S-sin(2*j(m)*pi/(n+1))/S)/((1+4*alpha*(1-cos(i(k)*pi/(n+1)))^2)*(1+4*alpha*(1-cos(j(m)*pi/(n+1)))^2));
    end
end
h = sum(((2*sin(j*pi/(n+1))/S - sin(2*j*pi/(n+1))/S).^2)./(1+alpha*(2-2*cos(j*pi/(n+1))).^2));
k1 = 2*alpha./(1-2*alpha*h); % the constant k1

% computation of matrix K2 with typical element (7)
i = 2:2:n;
j = 2:2:n;
K2 = zeros(n,n);
for k = 1:length(i)
    for m = 1:length(j)
        K2(i(k),j(m)) = (2*sin(i(k)*pi/(n+1))/S-sin(2*i(k)*pi/(n+1))/S)*(2*sin(j(m)*pi/(n+1))/S-sin(2*j(m)*pi/(n+1))/S)/((1+4*alpha*(1-cos(i(k)*pi/(n+1)))^2)*(1+4*alpha*(1-cos(j(m)*pi/(n+1)))^2));
    end
end
d = sum(((2*sin(j*pi/(n+1))/S - sin(2*j*pi/(n+1))/S).^2)./(1+alpha*(2-2*cos(j*pi/(n+1))).^2));
k2 = 2*alpha./(1-2*alpha*d); % the constant k2

% weights from Theorem 1
Finv_new = T*Lai*T + k1*T*K1*T + k2*T*K2*T;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Comparison with the weights from the HP filter based on the direct
% inversion of the matrix in (2)
p = n:-1:1;
P = eye(n);
P = P(p,:);
Q = diag(2*ones(n,1))+diag(-1*ones(n-1,1),1)+diag(-1*ones(n-1,1),-1);
g = [-2;1;zeros(n-2,1)];
A = eye(n)+alpha*Q*Q;
F = A-alpha*g*g'-alpha*P*g*g'*P;
Finv_old = inv(F); % Finv_old should be identical to Finv_new

